\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\newcommand{\matrx}[1]
{
	\ensuremath
	{
		\left (
		\begin{matrix}
			#1
		\end{matrix}
		\right)
	}
}

\begin{document}
	\section{Crank-Nicolson法与一维扩散方程}
	\footnote{参考：Sauer 《数值分析》。本笔记使用AI辅助。}
	我们回到耳熟能详的一维扩散方程：
	\begin{equation}
		u = u(x, t), 
		\left \{
		\begin{aligned}
			\pdv{u}{t} &= D \pdv[2]{u}{x}\qquad \text{扩散方程}\\
			u(x,t_0) &= f(x) \qquad \text{初始条件}\\
			u(x_0,t) &= 0 \qquad \text{左侧边界条件}\\
			u(x_n,t) &= 0 \qquad \text{右侧边界条件}\\
			x & \in [x_0, x_n]\\
			t & \ge t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	在隔壁笔记中，我们已经知道了扩散方程的显Euler迭代格式，其中假定$i=0,1,\cdots,n$共$n+1$个格子：
	\begin{equation}
		\text{显式} \quad \frac{u^{(k+1)}_i - u^{(k)}_i }{\Delta t} = D \frac{u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1}}{(\Delta x)^2}\\
	\end{equation}
	接下来引入今天的主题：隐式法和Crank-Nicolson法，
	Crank-Nicolson法的稳定性和数值精度比显式法更高（高稳定性是隐式法的经典特性）：
	\begin{equation} \label{eq_cn}
		\begin{aligned}
			\text{隐式} \quad & \frac{u^{(k+1)}_i - u^{(k)}_i}{\Delta t} = D \frac{u^{(k+1)}_{i+1} - 2u^{(k+1)}_{i} + u^{(k+1)}_{i-1}}{(\Delta x)^2}, \\
			\text{Crank-Nicolson} \quad & \frac{u^{(k+1)}_i - u^{(k)}_i}{\Delta t} = D \cdot \frac{1}{2} \left( \frac{u^{(k+1)}_{i+1} - 2u^{(k+1)}_{i} + u^{(k+1)}_{i-1}}{(\Delta x)^2} +  \frac{u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1}}{(\Delta x)^2} \right) \\
		\end{aligned}
	\end{equation}
	差别在于右侧空间导数项的时刻：
	\begin{itemize}
		\item 在显式法中，我们使用当前时刻$(k)$的场量估算空间导数值，这是非常常规的思路；
		\item 而在隐式法中，我们使用下一时刻$(k+1)$的场量估算空间导数值；
		\item Crank-Nicolson法则使用了这二种导数的平均。
		某种意义上说，Crank-Nicolson法类似于ODE的梯形Euler法。
		很奇怪，如此浅显的推广似乎直到1940s才被正式探讨
		\footnote{J. Crank, P. Nicolson, "A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type" Proc. Cambridge Philos. Soc. , 43 (1947)}。
	\end{itemize}
	但这里有一个非常显而易见的问题：我们如果都不知道下一时刻的场量，又如何计算它的空间导数？
	这是一个堪比先有鸡还是先有蛋的难题啊。
	在非线性ODE中，这往往需要迭代求解一个非线性方程；
	而幸运的是，扩散方程刚好是一个线性PDE，允许我们采取比较简单的做法。
	简单地整理\formula{eq_cn}，将$u$按时刻$(k)$分到两端（以Crank Nicolson为例），得到
	\begin{equation}\label{eq_cn2}
		u^{(k+1)}_i - \frac{1}{2}\alpha \left( u^{(k+1)}_{i+1} - 2u^{(k+1)}_{i} + u^{(k+1)}_{i-1} \right)
		= u^{(k)}_i + \frac{1}{2}\alpha \left( u^{(k)}_{i+1} - 2u^{(k)}_{i} + u^{(k)}_{i-1} \right)
	\end{equation}
	其中$\alpha = \frac{D \Delta t}{(\Delta x)^2}$。

	为了进一步说明，我们引入一点矩阵的语言。
	在隔壁Laplacian的笔记中，我们知道场量$\bvec u$可以写为向量的形式并对其运用$L$矩阵以求数值二阶导：
	\begin{equation}
		L \cdot \bvec u = 
		\begin{bmatrix}
			1 & 0 & 0 & 0 & \cdots & 0 \\
			1 & -2 & 1 & 0 & \ddots & 0 \\
			0 & 1 & -2 & 1 & \ddots & 0 \\
			\vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\
			0 & \cdots & 0 & 1 & -2 & 1 \\
			0 & \cdots & 0 & 0 & 0 & 1
		\end{bmatrix}
		\begin{bmatrix}
			u_0 \\ 
			u_1 \\ 
			u_2 \\
			\vdots \\ 
			u_{n-1} \\
			u_n
		\end{bmatrix}
		=
		\begin{bmatrix}
			u_0 \\ 
			u_0 - 2u_1 + u_2 \\ 
			u_1 - 2u_2 + u_3 \\
			\cdots \\ 
			u_{n-2} - 2u_{n-1} - u_n \\
			u_n
		\end{bmatrix}	
	\end{equation}
	（严格来说，相差一个$\frac{1}{(dx)^2}$）。
	这个矩阵的首尾行体现了固定边界条件，即保持$u_0, u_n$不变，
	若初始的$u_0=u_n=0$，即相当于固定边界为0。	
	于是，\formula{eq_cn2}现在也具有矩阵方程组的形式：
	\begin{equation} \label{eq_cn3}
		(I - \frac{\alpha}{2} L ) \bvec u^{(k+1)} = (I + \frac{\alpha}{2} L ) \bvec u^{(k)}
	\end{equation}
	之所以能这么做，是因为之前提及的“扩散方程是线性的”。
	其中$I=diag(1,1,...)$是单位阵。压缩\formula{eq_cn3}为更精炼的形式：
	\begin{equation}\label{eq_cn4}
		A \bvec u^{(k+1)} = B \bvec u^{(k)} 
		\qquad
		\text{其中}
		\qquad 
		A = I - \alpha/2 L \qquad	B = I + \alpha/2 L 
	\end{equation}
	即
	$$
	\bvec u^{(k+1)} = A^{-1}B \bvec u^{(k)}
	$$
	耶！
	
	然而，$L, A, B$和都是一个很大的矩阵$((n+1)\times(n+1))$，但只有对角线附近大约共$3n$个元素不为零，
	因此$L, A, B$等都是典型的稀疏矩阵。
	直接计算$A^{-1}B$会导致消耗大量内存并耗费较长计算时间。
	聪明的做法是，将$A, B$保存为稀疏矩阵，并使用稀疏乘法和稀疏方程组求解器来求解 \formula{eq_cn4}。
	并且，若$n$很大，可以考虑使用gmres这类近似算法。
	Python代码示例：
	\begin{center}
		u1 = scipy.sparse.linalg.gmres(A,B.dot(u0))[0] \#A,B是稀疏阵，例如scipy.sparse.csr\_matrix
	\end{center}
	(或许还可以尝试三对角阵的追赶法，可能更快)因此，更准确地说，Crank-Nicolson法是一个循环求解稀疏矩阵方程组的过程。
	
	我们容易将$A,B$推广以改变不同时刻下数值导数的权重，或者说，含量：
	$$
	A = I - \alpha \theta L \qquad	B = I + \alpha (1-\theta) L 
	$$
	其中$0 \le \theta \le 1$。当$\theta = 0$时，相当于显式法，当$\theta = 1/2$时是Crank-Nicolson法，当$\theta = 1$时是隐式法等。
	
	只要将$L$推广至二维或者三维形式，即可使用Crank-Nicolson法求解二维或者三维扩散问题。参考隔壁Laplacian方程笔记。
	
	Crank-Nicolson法还可以用于求解Schrodinger方程等。

\end{document}
